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ABSTRACT 

We describe a morphological imprint of magnetization found when considering the relative orien- 
tation of the magnetic field direction with respect to the density structures in simulated turbulent 
molecular clouds. This imprint was found using the Histogram of Relative Orientation (HRO): a 
new technique that utilizes the gradient to characterize the directionality of the density and column 
density structures on multiple scales. We present results of the HRO analysis applied to three models 
of molecular clouds in which the initial magnetic field strength is varied, but an identical initial tur- 
bulent velocity field is introduced, which subsequently decays. The HRO analysis was applied to the 
simulated data cubes and mock-observations of the simulations produced by integrating the data cube 
along particular lines of sight. In the 3D analysis we describe the relative orientation of the magnetic 
field B with respect to the density structures, showing that: (1) The magnetic field shows a preferen- 
tial orientation parallel to most of the density structures in the three simulations. (2) The orientation 
of B changes from parallel to perpendicular to the density structures in regions with density over 
a critical density tit in the highest magnetization case. (3) The change of relative orientation with 
respect to density is largest for the highest magnetization and decreases in lower magnetization cases. 
This effect is also present in the projected maps, where we consider the relative orientation of column 
density structures and polarization pseudovectors. The described change in the relative orientation of 
the magnetic field is potentially measurable using submillimeter polarization observations of molecular 
clouds and can be used to complement the estimates of magnetic field strength found with techniques 
based on the Chandrasekhar-Fermi method. In conjunction with simulations HROs can be used to 
establish a link between the observed morphology in polarization maps and the physics of molecular 
clouds. 

Subject headings: ISM: magnetic fields, ISM: clouds, polarization, submillimeter 



1. INTRODUCTION 

The study of the magnetic field in molecular clouds is 
crucial to understand the physical processes relevant in 
mole c ular cloud evolution and star formation (|Shu et al.l 
[19871 : iMcKee fc O striker! [20071 : lCrutcherll2012l ). A long- 
standing technique used to infer the morphology of the 
magnetic field projected on the plane of the sky is the 
measurement of polari zation of visible and near-infrared 
light background stars (iHiltner 1949; D avis fc Greensteinl 
1 19511 : iHeiles fe Haverkornll2012[ ). This technique assumes 
that aspherical dust grains are aligned with respect to 
the interstellar magnetic field and the observed position 
angle of linear polarization is parallel to the plane-of-t he- 
sky projection of the field (|Davis fc Greensteinl 1 19511 ). 

A more recent complementary technique is the mea- 
surement of the thermal polari zed emission of dust in 
the submillimeter wavelengths rtHildebra nd et all 119841 : 
iNovak et all 119971 : IVaillancourtl I2007D . Measurable de- 
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grees of polarization on the order of a few percent 
are typical for interste llar clouds in the submillime- 
ter a nd the far-infrared (jBenoit et al.ll2004l : IVaillancourtl 
120111 ). The mechanism which produces the alignment 
of the grains is still th e subject of active research 
(jLazarian fc Hoangl |2007| ). However, there is firm evi- 
dence indicating that it results from radiative torques 
produced by anis otropic radiation flux wi th respect to 
the magnetic field ([Hoang fc Lazari an 2008). Theoretical 
understanding of the grain alignment indicates that even 
in relatively dense clouds (visual extinction Ay ~ 10), 
dust polarization t races the underlying magnetic field 
(jBethell et al.ll2007D . 

The strength of the magnetic field has been esti- 
mated from observations of the projected component 
on the pl ane of the sky using the Chandr a sekhar-Fermi 
meth od ( Chandrasekhar fc Fermil 119531 : iHoude et al.l 
12001 INovak et al.l 12009) and comp aring the observe d 
morphology with theoretical models (K och et al.ll2012"l ). 
In parallel, significant efforts have been made to re- 
late the magnetic field morphology inferred from po- 
larization tothe obs^rved_s_tri^ both observed 
(|Li et alJjgOOl [C hapma n et al.l l20lTh and simu lated 
clouds ([Ostriker et al . 2001; Nakamur afc Lill2008l ). In 
the case of pre-stellar cores, a number of observations 
reveal an hourglass morphology in th e magnetic field 
(jGirart et al.l 120061 : lAttard et all 12009). and in the case 
of larger structures, the morphological analysis relies on 
visual inspection and st atistics of the orientation angles 
([Goldsmith et al.l l2QQ8h . The results of observations 
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are inconclusive in establishing a global trend of rela- 
tive alignment between the magnetic field and the struc- 
ture (|Goodman et al.lll99Ql ). However, theory and simu- 
lations establish very clear scenarios in which the relative 
orientation of the field reveals the role of magnetism in 
the formation of structures the cloud (|Fiege fc Pudritzl 
[200(1 [Nakamura fc Lill2008l ). 

Recent observations by the ESA Herschel Mis- 
sion reveal ub iquitous filam e ntary s tructures in molec- 
ular clouds (iPeretto e _jaL| 120121 : lAndre et all 120101 : 
iMolinari et al.ll201Q[ ). These observations together with 
starlight polarization measurements have been used to 
investigate the physical conditions of gravitational col- 
lapse and star formation within the filaments. Seminal 
examples in which the inferred magnetic field direction in 
adjacent lower density material i s perpendicular to the 



densest filaments are in Taurus (Chapman e t al. 
IPalmeirim et al.l 120131 ), Serpens (Sugita ni et al.l _ 
and the Pipe Nebula (jAlves et al.l I2008D . These stud- 
ies provide insight into the physical processes involved 
in the formation of filamentary structure but they do 
not probe the complex relationship of column density 
and magnetic field morphologies in the dense regions of 
molecular clouds, where stars predominantly form and 
which are fundamental to e stablishing a glob al picture of 
molecular cloud formation ([Lada et al.ll201Q[ ). 

The advent of new instruments for measuring 
the thermal dust po larized emissi on such as ALMA 
(jPeck fc Beas lev 2008 s ). BLASTPol ftP^^rt^ M^ 
and Planck ([Planck Collaboration et al.l 1201 if ) which 
measure polarization in scales ranging from pre-stellar 
cores to entire molecular cloud complexes, motivates a 
common scheme for the morphological analysis of polar- 
ization maps. It is in this context that we introduce the 
Histogram of Relative Orientation (HRO), a statistical 
tool for characterizing column density morphologies on 
multiple scales. This technique, developed for pattern 
recognition in machine vision, uses the gradient to char- 
acterize the directionality of the density structure in a 
2D map and here in a 3D simulated cube. It provides 
a robust characterization of the density field that can 
be used with polarization observations to investigate the 
relative orientation of density structure with respect to 
the magnetic field. 

This paper is organized as follows: Section [2] introduces 
the method of Histogram of Relative Orientation and 
implementation details of the technique are described in 
Appendix [A] Section [3] introduces the simulated molec- 
ular clouds used to characterize our statistical tool and 
details are given in Appendix [B] Section 2] presents the 
HRO analysis of 3D simulated clouds and Section [5] shows 
how the proposed diagnostic obtained in 3D can also be 
found in the projected 2D maps. Finally Section [6] sum- 
marizes the results and discusses the application of HROs 
to observations of real molecular clouds. 

2. THE HISTOGRAM OF RELATIVE 
ORIENTATIONS 

The Histogram of Relative Orientation (HRO) is 
inspired by a family of detection algorithms usu- 
ally called Histogram s of Oriented Gradients (HOGs, 
iLeonardis et "all [20061) . HOGs are used in pattern recog- 
nition and are developed in the context of machine vision. 
These algorithms use the gradient to describe the orien- 



tation of the edges in an image. The gradient is based 
on the difference between pixels in a neighborhood and 
quantifies the edges: the magnitude carries the informa- 
tion on the difference (edge strength) and the orienta- 
tion defines the direction perpendicular to a contour line 
(edge direction). A histogram of the orientation angles of 
the gradient rotated by 90° therefore summarizes statisti- 
cally the orientation of the features in an image. Further 
refinement of this method includes weighting the contri- 
bution of each pixel (or voxel in 3D) to the histogram 
by the magnitude of the gradient at each point, allow- 
ing characterization of the image using only the stronger 
edges. 

In HROs, the computation of the gradient is performed 
using Gaussian derivative kernels. The size of the Gaus- 
sian determines the area of the vicinity over which the 
gradient will be calculated. Varying the size of the Gaus- 
sian kernel enables the sampling of different scales and 
reduces the effect of noise in the pixels. 

The HROs use the gradient of the density in 3D simu- 
lated cubes to generate a vector field which characterizes 
isodensity contours and can be directly compared to the 
magnetic field vectors. In 2D, the gradient of the column 
density (VE) is compared to the orientation of the lin- 
ear polarization field, which is a pseudovector since it is 
invariant to a ±180° rotation. 

The implementation of the HROs is discussed in detail 
in Appendix [XJ 

3. MODEL PARAMETERS 

The cloud models used to characterize HRO are 
the result of magneto-hydrodynamical s imulations pro- 
duced with the RA MSES-MHD code ([Tevssierl 120021 : 
iFromang et al.l [20061 ). These clouds are isothermal with 
T = 11.44 K (sound speed c s =0.2 km s _1 ) for regions 
with density less than 10 4 times the initial density no 
and adiabatic for larger densities. The solutions are ob- 
tained in a cubic box of side L with initial uniform den- 
sity no = tih 2 = 536.41 cm -3 , which is comparable to 
the density in the Tau rus- Auriga Dark Cloud Complex 
(jStahler fc Pallal 120051 ). The boundary conditions of the 
box are periodic. 

The relative importance of gravity and thermal pres- 
sure forces is related to L and no- For this analy- 
sis we have chosen L = 4 pc; the Jeans length Lj = 
c s (tt /Gpo) 1 / 2 is 0.88 pc for po = A^o and \i = 2.4 m p . 
The total mass in the simulated cubes is 2.04 x 10 3 M 
and with L = 4 pc, it is comparable to small re- 
gions within Da rk Cloud Complexes such as Lupus I 
(|Cambresvl[T999h . 

Turbulence in the simulation is introduced as an ini- 
tial isotropic random velocity field 5v with a Kolmogorov 
spectrum (see Appendix [B]). This velocity field is not 
driven and decays during the cloud evolution. All 
simulations are initiated with kinetic energy Ek = 
100 pL 3 c 2 corresponding to an initial Mach number M. = 
cr v /c s =10 and include gravity. 

Three different initial magnetization cases are consid- 
ered and parameterized by the ratio of the isothermal 
sound speed to the initial Alfven speed or equivalently, 
the ratio of the plasma pressure to the initial magnetic 
pressure: /3 = c 2 /v\ = pc 2 s / (Bq / 'An) . We consider a 
"quasi-hydrodynamical" model with f3 = 100.0, an "in- 
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TABLE 1 

Model parameters and times of 
snapshots 



p 


< B > 


Snapshot Age 




( M G) 


(Myr) 


(V) (*/ b ) 


100.0 


0.35 


0.48 


0.11 0.25 


100.0 


0.35 


1.16 


0.27 0.60 


1.0 


3.47 


0.69 


0.16 0.36 


1.0 


3.47 


1.32 


0.31 0.68 


0.1 


10.97 


0.61 


0.14 0.32 


0.1 


10.97 


1.15 


0.27 0.59 



a Characteristic gravitational time scale tg — 
4.27 Myr 

k Flow crossing time scale tf — 1.95 Myr 



termediate magnetization" model with j3 = 1.0, and a 
"high magnetization" model with f3 = 0.1. The physical 
value of the magnetic field is given by: 




For the chosen initial density and temperature, the 
corresponding uniform magnetic field strengths are 0.35, 
3.47, and 10.97 fiG. The evolved fields are spatially 
nonuniform and can differ greatly from the initial values 
although the volume-averaged magnetic field is a con- 
stant Box in time. The value of f3 is proportional to the 
the mass-to-magnetic flux ratio in the simulation and 
cannot change with time. 

1/2 

The initial Alfven Mach number Ma = ( v<2 / v a) 
= Mp l/2 is 100, 10, and 3.16 for /? = 100.0, 1.0, 
and 0.1 respectively, and therefore all models consid- 
ered have initially supersonic and super-Alfvenic flow. 
The initially uniform clouds are unstable to compres- 
sions transverse to the mean magnetic field ("magneto- 
Jeans unstable") when the magneto-sonic wave cross- 
ing time excee ds the characteristic gravitatio nal contrac- 
tion time, t 9 ( Chan drasekhar fc Fermi|[l953[ ). All mod- 
els would be unstable by the magneto- Jeans criterion 
(L/Lj = 4.56 < /3 -1 / 2 ). The three simulations are ini- 
tially supercritical, with the mass M to magnetic flux <$> g 
ratio under the critical value [M /$ g ) C rit — ^/i^G 1 / 2 ), 
and (M/$g)/(M/$ g ) cr it « t g v A /{nL) is equal to 0.01, 
0.07, and 0.22. 

In Table [TJ we describe the simulation snapshots under 
consideration in terms of the sound crossing time t v = 
L/c s = 19.5 Myr, which is fixed owing to the isother- 
mal equation of state, the flow crossing time tf = L/<j v 
= 1.95 Myr, where we use the Mach number associated 
with the initial turbulent velocity and the characteris- 
tic gravitational time scale, t g = (ir/Gp) 1 / 2 = 4.27 Myr. 
Two snapshots in this study are taken: one at t\ ~ 0.03 t v 
and another at t^ y 0.06 t v . As in p receding studies 
(jOstriker et al.l200HlHeitsch et al.l2001l ). this paper con- 
centrates on structures that form as a consequence of tur- 
bulence, before gravity becomes important. The snap- 
shots are taken at t\ < £2 < tg when shock fronts moving 
through the gas initiate the formation of filaments and 
knots. 

Further details of the numerical methods used in the 
simulations are given in Appendix [Bj 
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Fig. 1. — Histogram of Relative Orientation showing the cosine of 
the angle between the magnetic field vector B and the gradient of 
the density Vn for the low, intermediate, and high-magnetization 
simulation cubes (/3 = 100, 1.0, and 0.1) in a snapshot taken at 
t ~ 0.06 t v . The histogram is normalized such that a random 
distribution of B and Vn would equal unity in each bin (black 
dotted line). The histograms calculated from the simulated cubes 
show a peak at cos (/> ~ which corresponds to the magnetic fields 
predominantly tracing the isodensity contours. 

4. HRO APPLIED TO SIMULATION CUBES 

The simulated cube contains the density, velocity, and 
magnetic field values for every voxel. The HROs are 
used to summarize the relative orientation of the mag- 
netic field vector B with respect to the gradient of the 
density Vn in each one of these voxels following the pro- 
cedure described in Section^ The histogram of the angle 
4> between two sets of random vectors in 3D is not uni- 
form, i.e., in 3D, two random vectors are more likely to 
be perpendicular than parallel. In 3D the uniformly dis- 
tributed quantity is cos</> and therefore we choose that 
quantity for the histogram. 

4.1. Results of the HROs 3D Analysis 

Figure [T] shows the HRO corresponding to all of the 
voxels in each of the three simulations for a snapshot at 
t2 ~ 0.3 t g . The peak at cos(/> = shows that Vn is 
predominantly oriented perpendicular to B in the three 
simulations. This corresponds to B mostly parallel to 
the isodensity contours in the simulated volume. 

The two snapshots considered here correspond to t < 
t g and t < t v . At these time scales, the dynamics of 
the gas are dominated by shocks and the over-densities 
are not formed by self-gravitation. The relative orienta- 
tion revealed by the HRO is the result of the magnetic 
field becoming strongly bent and stretched. The imprint 
of the initial magnetization was found by studying and 
comparing the relative orientation in the highest densi- 
ties regions. In these regions, the relative orientation is 
the consequence of the locally strong self-gravity and the 
interaction of the gas and the magnetic field. 

The colored curves in Figure [2] correspond to the HRO 
of regions in a particular density bin. The red and ma- 
genta curves correspond to regions with the highest den- 
sities, the gray and blue curves to regions with densities 
close to the mean density of the cube (logn = 2.73) and 
all others to intermediate densities. The algorithm used 
to seg ment t he cube in these density bins is described in 
Section IA.3I 
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Fig. 2. — HROs corresponding to simulated cubes with /3=100.0 (top), 1.0 (middle), and 0.1 (bottom) in snapshots taken at t ~ 0.03 t v 
(left) and t ~ 0.06 £ v (right). The colored curves within each plot correspond to voxels in the density ranges indicated in the figure. The 
histograms in the low magnetization case (top) peak at cos cf) ~ in regions with densities n > n which corresponds to the magnetic field 
(B) predominantly tracing the isodensity contours even at the greatest densities. The histograms from the intermediate (middle) and high 
magnetization (bottom) cases also peak at cos cf) ~ in regions with densities n ~ n. However, the histograms flatten in higher density 
regions. At the highest densities the histograms peak at cos</> ~ ±1. This corresponds to B tracing the isodensity contours in regions with 
n ~ n, then showing no particular relative orientation in intermediate density regions, and being oriented perpendicular to the isodensity 
contours in the highest density regions. 



Figure [2] illustrates the difference in the relative ori- 
entation of Vn and B in low and high density regions 
in simulations with low and high magnetization. In the 
low and intermediate magnetization cases (/? = 100 and 
P = 1 ), Vn and B are preferentially perpendicular to 
each other, which corresponds to B parallel to the iso- 
density contours. In the case of greater magnetization 
(/? = 0.1), Vn is mostly perpendicular to B at densi- 
ties close to the mean but the relative orientation pro- 
gressively changes when considering regions with greater 
densities until Vn and B are predominantly parallel to 



each other. 

With this progressive change in relative orientation, 
the HRO curve changes from convex to concave. To 
quantify this behavior we define the histogram shape pa- 
rameter: 

C = A c - A e (2) 

A c is the area under the central region of the HRO 
curve (—0.25 < cos</> < 0.25) and A e is the area in the 
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extremes of the HRO curve (—1.00 < cos</> < —0.75 and 
0.75 < coscj) < 1.00). This parameter characterizes a 
curve peaking at cos 0^0.0 (convex) as ( > 0.0 whereas 
a curve peaking at cos</> ~ ±1.0 (concave) corresponds 
to £ < 0.0 and a flat distribution corresponds to £ ~ 0. 
The uncertainty in the determination of £ is related to 
the standard deviation around the calculated area in each 
region, i.e., 

= yl°A c +*A e ( 3 ) 

With these definitions, a noisy HRO will either pro- 
duce ( ~ which is consistent with a random relative 
orientation or will increase cfq to a level in which it de- 
scribes no preferential orientation. 

Figure [3] shows the value of £ for the different column 
density segments in all six simulations. The red curve 
illustrates the progressive change of the relative orienta- 
tion of Vn and B with density in the simulations with the 
greatest magnetization: it changes from B preferentially 
parallel to the isodensity contours (£ > 0) for densities 
close to the mean (logn [cm -3 ] ~ 2.7) to B perpendicu- 
lar to the isodensity contours (£ < 0) in the regions with 
densities higher than about 6n. The curves correspond- 
ing to low and intermediate magnetization also show a 
similar change in the relative orientation of Vn and B 
with density although it is less pronounced. Thus we find 
that the slope of ( as a function of the mean density in 
the density bin is a parameter which characterizes the 
initial magnetization state of the simulated cloud. 

Additionally, we find that £ becomes negative (B is 
preferentially perpendicular to the isodensity contours) 
at the highest densities in the highest magnetization case. 
The value of the threshold density, ny, above which Vn 
and B are predominantly parallel (£ becomes negative) 
is also a parameter which characterizes the magnetiza- 
tion state of the simulated cloud: higher magnetization 
corresponds to lower values of tit- 

4.2. Discussion 

We have found two features of the relative orientation 
of Vn and B which characterize the initial magnetiza- 
tion: 1. The rate of change from perpendicular Vn and 
B to parallel in regions with increasing mean bin den- 
sity. 2. The value of the density nr over which Vn and 
B are parallel. We observe that the rate of change of £ is 
greater and that nx is lower with higher magnetization. 
These effects are related to the balance between the mag- 
netic forces, the turbulence and the gravitational forces 
in each simulation. 

The three simulations are initially supercritical. Su- 
percritical clouds can collapse both parallel and perpen- 
dicular to the field with unlimited asymptotic density, as 
opposed to subcritical clouds which can collapse along 
the field but not perpendicular to it, forming flattened 
structures where the peak density would be limited by 
the thermal pressure. 

The supersonic perturbations resulting from the initial 
turbulence are amplified or suppressed depending on the 
strength of the initial magnetic field (Draine & McKed 
|1993[ ). The dynamical importance of the magnetic field 
before self-gravity becomes important can be estimated 
by comparing the initial magnetic energy density B 2 /8tt 



and the initial kinetic energy density pa 2 / '2. This ratio is 
(v A /cr v ) 2 = l/(M 2 /3) and it is equal to 0.001, 0.1, and 1.0 
for the low, intermediate, and high magnetization. In the 
low magnetization case, the Alfven speed va < c s and the 
kinetic energy density is dominant. In consequence, the 
fields become strongly bent by motions of the matter. 
In the high magnetization case, the magnetic and the 
kinetic energy are closer initially, and with va > c s , the 
matter is constrained to move along stiff magnetic field 
lines. 

The relative orientation between supersonic shock 
fronts and the magnetic field could in principle be the 
mechanism which produces over-densities with a partic- 
ular orientation with respect to B. Alignment between 
B and the gas velocity field h as been reported in previous 
works (|Baneriee et alJ[2Q09h : however the exact process 
which causes this relative orientation is not well under- 
stood yet. The relative orientation of B perpendicular 
to the over-dens i ties is close to the scenario described 
by iMouschoviasI (|1976l ). However, the interaction be- 
tween supersonic turbulence and magnetic field which 
produces the final configuration of B and density struc- 
tures is still the subject of active research ( Klessen 20111 : 
iSchneider etal]l2Ql"lt iHennebelle fc Fakaronell2012[ ). 

In the low-magnetization case, Vn is perpendicular to 
B (B along isodensity contours) even in the regions with 
the highest densities. In this case, the gravitational forces 
dominate over the magnetic field ((M/ <& g )/ [Mj & g ) C rit ~ 
0.01) and in the densest regions the decrease of £ is re- 
lated to the isotropic gravitational collapse breaking the 
alignment of Vn perpendicular to B, leading to ( ~ in 
the densest regions. 

In the high-magnetization cases, Vn is parallel to B 
(B perpendicular to the isodensity contours) in the re- 
gions with the highest densities. In this case, the gravi- 
tational forces are less dominant ((M '/ g ) / \M / g ) cr it ~ 
0.22) and the channeling of the flow along the magnetic 
field is reinforced by the perpendicular magnetic sup- 
port against gravitational collapse, leading to flattened 
structures perpendicular to the magnetic field which cor- 
respond to ( < in the densest regions. 

5. HRO APPLIED TO OBSERVATIONS OF THE 
SIMULATION CUBES 

The HRO analysis is applied to the maps of the col- 
umn density E and polarization obtained by integrating 
the six simulated data cubes along three different lines of 
sight (x, y, or z-axis). Details of the mode l of projection 
using Equation IB II are given in Appendix IB. 11 The pro- 
jection of the Stokes parameters into the plane of the sky 
incorporate s a g rain alignment efficiency model described 
in Equation IB2I This uses the density in each voxel as a 
proxy for what might actually be a column density de- 
pendence of the grain alignment mechanism. Here we 
compare results for uniform alignment efficiency (e = 1) 
with the results for a model with e = 1 at n < no but 
decreasing e with increasing n in regions with n > no- 

The grain alignment is locally such that the polariza- 
tion vector is perpendicular to the magnetic field. The 
projected polarization is therefore perpendicular to a 
weighted projection of the magnetic field, a pseudovector 
that we call Bpos 1 - 

The HROs of the simulated maps are calculated using 
a process analogous to that used in the 3D data cubes. 
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Fig. 3. — HRO shape parameter £ which parameterizes the relative orientation of the magnetic field (B) and the gradient of the density 
Vn. C > 0-0 corresponds to a HRO showing B predominantly perpendicular to Vn (B parallel to the isodensity contours). £ ~ 0.0 
corresponds to a flat HRO showing no predominant relative orientation between B and Vn. £ < 0.0 corresponds to a HRO showing B 
predominantly parallel to Vn (B perpendicular to the isodensity contours). The HROs of the low magnetization case show B predominantly 
parallel to the isodensity contours with £ > 0.0 even in the higher density regions. In contrast, the HROs of the high magnetization cases 
show B parallel to the isodensity contours with £ > 0.0 in the low density regions, changing to B perpendicular to the isodensity contours, 
with £ < 0.0 in the highest density regions. The shape parameter £ as a function of density decreases faster in higher magnetization 
simulations. 



The histogram of the angle <\> between two sets of randon 
vectors in 2D is uniform, and therefore we choose (j) as 
the variable. 

5.1. Results of the HROs 2D Analysis 

Figure 2] shows the HROs corresponding to the an- 
gle (j) between the projected magnetic field pseudovector 
Bpos and the gradient of the column density VE in the 
projected maps. The peak at 90° reveals that the rela- 
tive orientation between Vn and B observed in the 3D 
(Figure [T]) is also present in the polarization and column 
density maps and Rpos predominantly follows the iso-E 
contours. This result can be confirmed by visual inspec- 
tion of the ~Bpos and E maps as shown in Figure 

As in Section HJ we divided the maps in column den- 
sity bins to check how the relative orientation changes 
in the highest density regions. Figure [6] shows the 
HRO curves corresponding to different column density 
regimes. The behavior of the HROs in 2D is very simi- 
lar to that observed in the 3D analysis: Rpos and VE 
are mainly perpendicular in all the E bins of the low- 
magnetization map. The HROs of the intermediate and 
high-magnetization cases show ~Bpos and VE predomi- 
nantly perpendicular in regions of the map with E ~ E 
and parallel in the regions of the map with the greatest 
E. 

In the same way as in the 3D analysis, the change 
in the HRO curves is parameterized by the histogram 
shape parameter £ defined in Equation 14.11 In this case 



A c is the area under the central region of the HRO curve 
(67.5° < (j) < 112.5°) and A e is the area in the extremes 
of the HRO curve (0° < <j> < 22.5° and 157.5° < < 
180°). Figure [3 shows the value of ( as a function of the 
central column density of each bin. The HROs of the low 
magnetization case show Rpos predominantly parallel to 
the isodensity contours with ( > 0.0 even in the highest 
E regions. In contrast, the HROs of the intermediate 
and high magnetization cases show ~Bpos parallel to the 
isodensity contours with ( > 0.0 but changing into ~Bpos 
perpendicular to the isodensity contours with ( < 0.0 
predominately in the highest E regions. 

5.2. Discussion 

There are two main effects to be considered in the re- 
covery of the imprint of magnetization on the relative 
orientation in the projected maps. First is the effect of 
the integration along the line of sight. The second is the 
effect of the alignment efficiency in weighting the recon- 
struction of the projected magnetic field. 

Figure [3 shows that the relative orientation in the pro- 
jected map is analogous to the one measured in 3D: 
Bpos 1 and VE are preferentially aligned perpendicular 
to each other. In 3D the highest density regions where 
we observed B perpendicular to the isodensity contours 
are surrounded by lower density shells where B is paral- 
lel to the isodensity. In 2D we observed Rpos parallel 
to VE in the highest density regions of the map show- 
ing that the relative orientation in the lower density shell 
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between the projected magnetic field pseudovector Bpog and the 
gradient of the column density VE in projections of the low, inter- 
mediate, and high-magnetization simulation cubes ((3 = 100, 1.0, 
and 0.1) in a snapshot taken a t t ~ 0.06 t v and using grain align- 
ment efficiency from Equation IB2I with p = (top) and p = — 1 
(bottom). The projections are obtained following Equation IB1I 
The histogram is normalized such that a random distribution of 
Bpo^ and VE would equal unity in each bin (black dotted line). 
The histograms calculated from the simulated cubes show a peak 
at cf) ~ 90° which corresponds to the magnetic fields predominantly 
tracing iso-E contours. 



does not dominate in the projected map. More detailed 
modeling of the integration along the line of sight and its 
effect on the HROs and ( will be the subject of future 
studies. 

An interesting issue for the line of sight integration is 
the relative orientation of the mean magnetic field with 
respect to the line of sight. All of the results presented 
here correspond to lines of sight perpendicular to the 
plane of the initial magnetic field (Y or Z for initial B 
along the X axis). Small variations about this orientation 
of the mean magnetic field with respect to the line of 
sight do not affect the conclusions of the HRO study. 
However, the HROs corresponding to projections along 
the X-axis (line of sight parallel to the initial magnetic 
field direction) are not distinctively different for different 
initial magnetization. Thus, further studies are required 
to draw conclusions from observations of clouds where 
the mean magnetic field is oriented very close to the line 
of sight. 

The effect of the alignment efficiency is illustrated in 
the two plots in Figure [7J The slopes of the £ — E 
curves corresponding to different initial magnetization 
are clearly distinguishable in the case of uniform e = 1. 
However, when considering a decreasing alignment effi- 
ciency, e oc n _1 , the curves corresponding to the low 
and intermediate magnetization are degenerate and only 



the curve corresponding to the highest magnetization is 
distinctive. 

We find that the slope of the ( — E curve is strongly 
dependent on the grain alignment mechanism, which is 
yet to be understood. In spite of that, £ < corresponds 
to magnetic field dominance. The column density Et 
at which the relative orientation changes can be used 
potentially to estimate the magnitude of ~Bpos- 

6. SUMMARY AND DISCUSSION 

We have introduced the Histogram of Relative Orien- 
tations (HRO), a tool to study the relative orientation of 
the magnetic field and the structure in MHD simulations 
and in polarization observations. The HRO quantifies 
the relative orientation of the magnetic field B and the 
isodensity contours. The contours are characterized by 
the gradient of the density Vn in the simulated volume 
or the gradient of the column density VE in projections 
of the simulations. 

Using HROs on a set of simulated molecular clouds 
with decaying supersonic turbulence and with low, inter- 
mediate, and high initial magnetization (/? = 100.0, 1.0, 
and 0.1) we found a significant imprint of the magne- 
tization in the relative orientation of the magnetic field 
and the density structures. We found that B is oriented 
predominantly parallel to the isodensity contours in the 
three simulations. However, on dividing the simulated 
region into density bins and considering the relative ori- 
entation in each bin, in the high magnetization simula- 
tion we found that at the highest densities B becomes 
perpendicular to the isodensity contours. 

We parameterize the change in relative orientation 
with the histogram shape parameter £. We show the 
behavior of £ in regions with different densities in the 
simulated volume (Figure [3]) and different column densi- 
ties in the projected map (Figure [7j). The study of this 
change in relative orientation can be applied to observa- 
tions of molecular clouds. Establishing a link between 
observations and MHD simulations will provide an op- 
portunity to study the role of the magnetic field in the 
formation of structures in molecular clouds. 

HROs can also be applied to observations where the 
polarization pseudovectors are only available on sparse 
regions of the map. Additionally HROs can use a range 
of gradient kernels to sample a particular range of scales 
and reduce the effect of the noise in the observed maps. 
With the aid of further simulations, the column density 
Ex at which the field changes its relative orientation may 
be used to estimate the strength of the magnetic field in 
the plane of the sky, complementing current estimates 
obtained with the Chandrasekhar-Fermi method. 

6.1. Relation to Existing Studies 

Studies of the morphology of the polarization obser- 
vations and its relation to the magnetic field strength 
in molecular clouds have gene rally been based on the 
Chandrasekhar-Fermi m ethod ([Chandrasekhar fc Fermil 
119531 : iHoude et al.ll2009f ). These studies relate the dis- 
persion of the polarization orientation to the magnetic 
field strength in the cloud, an effect whic h is evident in 
the previou s studies of MHD simulations (|Q striker et al.l 
|2001[ Il999f ) and in the polarization maps obtained from 
the simulation set in this paper (Figure [5]) . 
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Fig. 5. — Maps of the logarithm of the column density (logS) and overlaid magnetic field pseudovectors as determined by Equations IB1I 
and IB2I with p = — 1 and no = 500 cm -3 . These maps correspond to projections along the z-axis of simulations with (3 = 100 (top), 1.0 
(middle), and 0.1 (bottom) in snapshots taken at t ~ 0.03 t v (left) and t ~ 0.06 t v (right). 
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Fig. 6. — HROs corresponding to simulated cubes with /3=100.0 (top), 1.0 (middle), and 0.1 (bottom) in snapshots taken at t ~ 0.03 t v 
(left) and t ~ 0.06 t v (right). The colored curves within each plot correspond to pixels in the column density ranges indicated in the figure. 
The HROs from the low magnetization simulation (top) peak at cf) ~ 90 in regions with column densities E > E which corresponds to the 
projected magnetic field (Bpos) predominantly tracing the iso-E contours even at the greatest column densities. The histograms from 
the intermediate (middle) and high magnetization (bottom) cases also peak at ~ 90 in regions with column densities E ~ E. However, 
the histograms flatten when considering higher E regions and in the highest densities they peak at (j> ~ 0° or 180°. This corresponds to 
Bpos tracing the iso-E contours in regions with E ~ E, then showing no particular relative orientation in intermediate E regions, and 
predominately orienting perpendicular to the iso-E contours at the highest column density regions. As in the 3D analysis, the chan ge o f 
the relative orientation inferred from the histogram in different column density regimes is parameterized by Q as defined in Equation 14.11 
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Fig. 7. — Histogram of Relative Orientation shape parameter £ which parameterizes the relative orientation of the projected magnetic 
field (Bpos) and the gradient of the column density VS. £ > 0.0 corresponds to a HRO showing Bpos predominantly perpendicular 
to VS (Bpos parallel to the iso-E contours). £ ~ 0.0 corresponds to a flat HRO showing no predominant relative orientation between 
Bpos and VE. £ < 0.0 corresponds to a HRO snowing ~Bpos predominantly parallel to VS (Bpos perpendicular to the iso-S contours). 
As in the 3D case, the HROs of the low magnetization case show TSpos predominantly parallel to the isodensity contours with £ > 0.0 
even in the higher density regions. In contrast, the HROs of the intermediate and high magnetization cases show Rpos parallel to the 
isodensity contours with £ > 0.0 at low densities and changing into Rpos perpendicular to the isodensity contours, with £ < 0.0 in the 
highest density regions. The shape parameter £ as a function of density decreases faster for the higher magnetization simulations. 



An alternative approach is the study of the relative 
orientation of the magnetic field inferred from polar- 
ization observations and the orientatio n of the column 
densit y structure. In a seminal work iGoodman et al.l 
(j!99Q[ ) considered the alignment of starlight polarization 
and extinction maps in six different molecular clouds 
and found no particular orientation which globally de- 
scribed the observations. Similar studies (iTassis et al.l 
12009b iSugitani et all 1201 H iPalmeirim et al.l 2013) have 
considered a qualitative description to relate the mag- 
netic field observations to the formation of the filamen- 
tary structure. HROs represent a quantitative tool to 
evaluate the relative orientation of the magnetic field and 
create a general scheme to compare observations of dif- 
ferent clouds with the results of MHD simulations. 

The observed morphology of the submillimeter polar- 
i zed emission from the dense reg ions in molecular clouds 
(|Vaillancourt & Matthewsll2012f ) shows a diversity of rel- 
ative orientations which can potenti ally be correlated 
with the magnetic field in each region. iKoch et al.l (|2012h 
have already found a correlation between the iso-S con- 
tours and the inferred magnetic field directions which is 
used to fit a model of the magnetic field. The HRO pro- 
vides a common scheme to relate these observations to 
larger scale observations which cover entire cloud com- 
plexes. 

6.2. HROs and Future Work 

In this article we have focused our studies to the case 
of an isothermal molecular cloud with a particular be- 
havior of turbulence and magnetization in order to char- 
acterize our method. Nevertheless, the HRO analysis 
can be extended to simulations i ncluding more complex 
realizations of molecul ar clouds (|Nakamura & Lil [20081 : 
iHennebelle et al.l l2QQ8h . Characterization using HROs 
allow us to link a large volume of polarization observa- 
tions with the great diversity of physics which can be 
modeled in MHD simulations. 

The HROs of the simulations discussed in this paper 



evaluate the relative orientation in scales ranging from 
2.34 x 10~ 2 pc to 3.83 x 10" 1 pc (kernels of 3 x 3 to 
49 x 49 pixels), which correspond to angular scales of 
~10" to ~10' for a cloud located at a distance of 150 pc. 

The results of the HRO analysis show that the imprint 
of magnetization on the molecular cloud is present in at 
least two diagnostics of the relative orientation of the 
magnetic field: 

1. The change in the relative orientation of B and the 
column density structures as a function of column 
density. 

2. The value of the column density Y*t at which B 
shifts from parallel to perpendicular to the iso-S 
contours. 

These two parameters allow the systematic comparison of 
simulations and observations; however, characterization 
of the magnetization conditions using these diagnostics 
should be accompanied by careful modeling of the grain 
alignment efficiency. 

The study of multiple ages, environments and 
dominant physics in simulations is one of the pri- 
mary challenges facing the construction of a com- 
plete and coherent picture of the star formation 
process. The ad vent of results from expe r iment s 
such as Planck dPlanck Collaboration et al.l 120111). 
ALM A (jPeck fc B easlev 2008IT PI LOT (jBernard et al.l 
l2010l). BLASTPol dPascale et al.l I20TI ). and SOFIA 
(|Vaillancourt et al.ll2007t ) will produce an unprecedented 
volume of observations of molecular clouds. In the future, 
the use of HROs and complementary tools will facilitate 
the systematic characterization of synthetic observations, 
and identify the common elements in the observations 
and models that reveal the role of the magnetic field and 
turbulence in the formation of structure in the interstel- 
lar medium. 
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APPENDIX 
CONSTRUCTION OF THE HRO 
The implementation of the HROs is made with the procedure described in the following subsections 

Calculation of the Gradient 
The orientation of the isodensity contours are characterized by the gradient of density: 



^sJ.^U^HW, 5 (A1) 

The subindex I is related to the size of the area on which the gradient is calculated or as subsequently shown, the size 
of the Gaussian derivative kernel. To zeroth order, the components of the gradient can be calculated using forward 
differences of adjacent pixels: 

dn 

— ~n( x + l,y)- n(x, y) (A2) 
ox 

Derivatives are linear and shift- invariant, and so the gradient calculation can be done using convolution. Forward 
differences can be calculated using the 2x2 Roberts kernels, but such a small kernel is too sensitive to noise. The 
calculation of the gradient can be improved by using the central difference: 

dn n{x + 1, y) — n{x — 1, y) dn n(x, y + 1) — n(x, y — 1) 

~d~x = 2 ' djj = 2 (A3) 

The corresponding kernels for convolution calculation are the 3x3 Prewitt kernels: 

d i/- 101 \ d l/- 1 - 1 -^ 

s-5(:!!i!)'*"5(!!?) (A4) 

Noise reduction can be achieved by averaging over the vicinity of each pixel. For the convolution method this is 
implemented to first order using the 3x3 Sobel kernels. Convolving with these kernels is equivalent to smoothing the 
image over a 3 x 3 region and then calculating first derivatives, which is a particular case of the Gaussian Derivatives 
method: 

sMziiliM^ 1 (A5) 

Gaussian Derivatives 

Convolution and differentiation are commutative and associative; therefore the smoothing and derivative operators 
can be written as: 

i^ =I *l G (A6) 

G is the 2D- Gaussian kernel and I is a 2D image. This equation implies that filtering an image and subsequently 
calculating the gradient is equivalent to the convolution of th e image with a kernel that is the first derivative of a 
gaussian kernel. This operation is called a Gaussian derivative (|YoundH986[ ). 

For the present study, we use t he G aussian derivatives method and following Equation I A6I we obtained each compo- 
nent of the gradient in Equation lA.ll bv convolving slices of the density cube Nk = n^k and N{ = riijk with a kernel 
formed by the derivative of a I x I two-dimensional Gaussian G: 

dn\ __ d „ 
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Fig. 8. — Left: Filamentary structure in the simulated column density. Overlaid are the gradient vectors obtained with derivative kernels 
with 9x9 (yellow), 25x25 (red), and 49x49 (green) pixels. The squares on the upper right corner of the image show the sizes of these 
kernels. Right: Histogram of orientation angles of the iso-X contours (ip) calculated with each derivative kernel. The histograms show that 
the structure is predominantly oriented at ip from 135° to 150° but also shows a secondary structure at 45° (ip is measured counterclockwise 
from the top of the map, using the convention for position angles of polarization). 



and in the 2D projected maps 



OX 



dx J 
dy 



(A8) 



dy 



Figure [8] illustrates the characterization of a simulated image using the gradient of the column density VS. The 
vectors plotted over the image show the direction of the gradient obtained using Gaussian kernels of different sizes. 
The histogram shows the orientation angles of the structures in the image, measured counterclockwise with respect 
to the y-axis of the image (convention for polarization position angles). Each peak in the histogram corresponds to 
a dominant orientation of the features in the image. A histogram of a completely random map would be flat. The 
histogram of a straight rod would be a delta function centered on its orientation angle. The histogram of a blob would 
have a peak corresponding to the orientation of its semi- major axis. 

Calculation of the Angle 

The angle <p between the gradient and the magnetic field at each voxel is calculated using a combination of the scalar 
and vector product of vectors: 

B X Vn ^ (A9) 



arctan 



B • Vn 



The result of this calculation is a cube (a map in 2D) with values of <\> at each voxel (pixel in 2D). As explained in 
Section 3J in the study of the alignment in 3D it is more natural to construct the HROs as a function of cos (j) and not 
directly (j). The angle between the isodensity contours and the magnetic field is (j) =b 90°. The histogram of values of 
cos 4> (4> in 2D) weighted by the magnitude of the gradient at each voxel (pixel) is what we call Histogram of Relative 
Orientations. 

Figure [9] shows a histogram of cos (j) calculated on a simulated cube and (j) calculated on a projection of the cube 
using three derivative kernels. Using these three kernels we observe the same preferential relative orientation between 
B and Vn in the cube and ~Bpos and VS in the projection. 

Segmentation 

Individual regions of an image can be studied by dividing the images into bins of the parameters that define it and 
subsequently masking each bin. This process is known in computer science as image segmentation. Segmentation is 
a subject of active research in computer science and its implementation in 3D (2D) HROs is limited to dividing the 
cube (image) in density (column density) bins with equal number of voxels (pixels). 

The objective of the segmentation is the study and comparison of the relative orientation in density regimes where 
the dominant physical processes and the dust alignment efficiency e are different. The segmentation by column density 
of the maps resulting from the projection of the simulated cubes is illustrated in the left-hand side of Figure [101 Using 
the column density distribution, which is close to log-normal, we produce density bins with equal number of pixels. 
An equal number of gradient vectors in each bin guarantees comparable statistics for each density bin. 
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Fig. 9. — HROs obtained by applying three different derivative kernels to the (3 = 0.1 simulation cube (left) and to a projected map 
obtained by integrating the cube along the z-axis (right). The similarity of the curves within each plot shows that the predominant relative 
orientation is present in regions with scales ranging from 3 pixels to 49 voxels (pixels). 



The segmentation by density or column density is motivated by the change of physical conditions in the densest 
regions with respect to less dense regions. In the densest regions self-gravity becomes relevant and the kinetic energy 
density is lower than the magnetic energy density, in contrast with the less dense regions where the balance between 
magnetic and turbulent energy is different. The segmentation of maps by polarization percentage, line of sight velocity 
or polarization angle dispersion is a potentially useful way to study physical processes in observations of molecular 
clouds and will be addressed in future works. 



MOLECULAR CLOUD SIMULATIONS 

The model clouds are created by integrating the compressible ideal MHD equations using the RAMSES-MHD code 
(|Tevssierll2Q0llFromang et al.l|20()6). RAMSES-MHD is a N-body and MDH code w ith 3D Adaptive Mesh Refinement 
(AMR). The refinement criteria are based on density as described in lDib et~a71 (j2010f ). RAMSES-MHD uses constrained 
transport to guarantee that V • B = to machine accuracy at all time and uses the MUSCL-Hancock scheme, a finite 
volume method that combines good accuracy with fast execution. The solutions are obtained in a cubic box of side L 
with grids of 512 3 zones, integrated from the AMR cubes with effective resolution 2048 3 . We apply periodic boundary 
conditions in all models, simulating a portion of the interior of a molecular cloud. 

The energy equation used in this simulation set is the barotropic equation of state: for densities under the critical 
density parameter (n c = 5.364 x 10 7 cm -3 ) the gas is isothermal with sound speed c s = 0.2 km s _1 , which corresponds 
to a temperature of 11.44 K. In the absence of fully time-dependent radiative transfer this represents a reasonable first 
approximation for the gas at vol ume densities hig her than the mean (comprising most of the matter) for conditions 
appropriate to molecular clouds (|Scalo et al.lll998f ). When the density reaches n c , the gas becomes adiabatic. 

The gravitational potential is computed from the density using standard Fourier methods. The periodic boundary 
conditions exclude the k = component of the density. Therefore, the gravitational potential obeys the modified 
Poisson equation V 2 (/>g = 47rG(p — p), where p = M/L 3 is the mean density (mass divided by the volume of the box). 

Turbulence in the simulation is introduced as a random velocity field 5v. £v is a Gaussian random perturbation 
field with a power spectrum \5v^\ 2 oc /c -11 / 3 , subject to the constraint V • 5v = so that the initial velocity is non- 
compressive. This spectrum corresponds to the Kolmogorov e nergy spectru m, E(k) oc /c -5 / 3 , and it is comparable to 
the spectrum inferred for large-scale cold interstellar clouds (jLarson| [T981) and the spectrum that naturally follows 
from the evolution of incompressible turbulence. An identical realization of the initial velocity field is used for all of 
the models, so that the initial states of the simulations differ only in the strength of the mean magnetic field. The 
turbulence decays during the cloud evolution (it is not driven). 



Projecting 3D Data Cubes with Magnetic Field to 2D Maps with Polarization 

The observed intensity of the polarized emission results from the combine d effect of the dust grain a lignment efficiency 
e and integration along the line of sight. The radiative torques mechanism (jLazarian fc Hoangll2QQ7| ) indicates that the 
degree of alignment depends on the exposure of the grains to an anisotropic starlight background. This dependence 
implies that the radiative torques mechanism should be more efficient in regions of relatively low extinction. 

The maps of the Stokes parameters I, Q, and U are produced by integrating the density n and the magnetic field B 
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Fig. 10. — Example of the map segmentation process introduced in Section [A. 31 Left: Column density distribution in projections of three 
simulations. Overlaid colors correspond to density bins with equal number of pixels used in HRO analysis. Right: Column density maps 
of the projections of three simulations. Overlaid colors correspond to regions with column densities shown in the histograms. 



in the cube along the k-direction (either x, y or z) according to: 

k 

Qij * £ ^ n m(4%Bg - BgBg) 



(Bl) 



Uij oc 2e. 



ijk' "iJK-^ijk ijk 
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In this study we use the density as a proxy of the column density dependence of e. For each cell of the simulated 
data cube we evaluate: 

/ (n ijk y if n ijk > n , R . 

e ^-\1.0 ifn ijk <n ^ 

In this toy- model, p < accounts for depolarization effects (Matt hews et al J 120011 : iHoude et al.ll2002l ) in high column 
density regions, p = is used to test the line of sight integration independently of the environment dependence of e. 
In this study we have chosen no = 500 cm 3 , corresponding to a relatively diffuse region within a molecular cloud. We 
present results obtained with p = and p = — 1. These choices intend to reproduce some of the effects observed in 
polarization maps. A more detailed calculation of e, including modeling the grain alignment mechanisms is beyond 
the scope of this work. 
The polarization percentage P is calculated as: 




This quantity, of particular observational interest (M atthews fc Wilson! 120021 : [Martin 2007), is shown for reference in 
Figure [11] as a function of column density E in projections of three simulations for two models of e, with p = and 
p = — 1. The different distributions seen there and the systematic behavior, specially the decline of P with increasing 
E in the right panels, offer diagnostics relating parameters in MHD simulations to observations. 
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Fig. 11. — Polarization percentage P in units of maximum polarization percentage Pmax as a function of column density obtained from 
the projections of the 3 simulat ions with uniform p olarization efficiency e = 1.0. Pmax ~ 20% is determined empirically from submillimeter 
diffuse Galactic dust emission (Benoit et al. 2004). These different distributions and the systematic behavior, especially on the right, offer 
diagnostics relating parameters in MHD simulations to observations. 



